####################################################################################

#What You See and What You Get
#This file conducts the analysis for W1 survey
#Code to create Tables, Figures, and data cited in the main text and appendix 
#using Routputs created by running the previous three files (see documentation)


####################################################################################

rm(list=ls())
library(AER)
library(tidyverse)
library(readxl)
library(lubridate)
library(sandwich)
library(xtable)
library(corrgram)
library(forcats)
library(estimatr)
library(fastDummies)
library(randomizr)
library(ANOVAreplication)
library(miceadds)
library(clubSandwich)
library(texreg)
options(scipen=6)

source("Code/functions.R")

########################### Main Paper

#Main ITT and CACE effects
load("Routputs/out-W1-estimates-pooled.RData")

######################### Figure 1(a)

itts <- cbind(out.all$out.itt.lm$out.lmc 
              ,out.all$out.itt.lin$out.linc)  
itts_tp <- itts[1:4,c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp

the.labels.plot_tp <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                         "Attribution to\nLula/Dilma/PT", 
                         "Attribution to Paes/\nPezao/Temer/MDB") 

png(file = "Figures/fig-1-a.png",
    width = 8, height = 11, units = 'in', res = 600)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(4,12,2,1))
plot(xs,c(-.75,-(nrow(tp)+2)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, 
     font = c(2, 1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1), 
     cex.axis = 1.3) 
abline(v=0,lty=1,col=gray(0.7))
abline(h=c(ys)+0.5,lty=3,col=gray(0.7)) 
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT+Ctrls (via Fixed Effects)","ITT+Ctrls (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.015),xpd=NA)
dev.off()


# Figure 2
load("Routputs/out-W1-estimates-pooled-hetwait.RData")

allcwf <- out.waitenroll[[2]]
tp <- allcwf[c(1:6, 8), ]

the.labels.plot_tp <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                         "Attribution to\nLula/Dilma/PT", 
                         "Attribution to Paes/\nPezao/Temer/MDB",
                         "Lula Evaluation", "Dilma Evaluation", "Paes Evaluation")

png(file = "Figures/fig-2.png",
    width = 8, height = 11, units = 'in', res = 600)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(4,12,2,1))
plot(xs,c(-.75,-nrow(tp)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects" 
     ,main="Recent Lotteries Survey"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, font = 2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-4.5),col=gray(0.7))
#abline(h=c(-19.5),col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT+Ctrls (Short Wait)","ITT+Ctrls (Long Wait)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.015),xpd=NA)
dev.off()


##################### Online Appendix

################### Appendix C

#Balance Tests
#Objects
load("Routputs/out-W1-balancepooled.RData")
load("Routputs/out-W1-balancebylottery.RData")

#Tables
pol.bal <- xtable(bind_rows(pol.bal.stats$balage[c(1:4, 8)], 
                            pol.bal.stats$balnoage[c(1:4,8)]))

wald <- xtable(pol.bal.stats$wald)
pol.bal.stats$p.permutation #p-value for Table C.1

print(pol.bal,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c3.tex",sep=""))

print(wald,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c1.tex",sep=""))

####### Balance by lottery
#Edital 17
ed17 <- xtable(bind_rows(bal.stats$edital17.2016$balage[c(1:4, 8)], 
                         bal.stats$edital17.2016$balnoage[c(1:4,8)]))

wald17 <- xtable(bal.stats$edital17.2016$wald)
bal.stats$edital17.2016$p.permutation #permutation p-value for Table C.5 edital 17

print(ed17,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c7.tex",sep=""))

print(wald17,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c5a.tex",sep=""))


#Edital 20
ed20 <- xtable(bind_rows(bal.stats$edital20.2016$balage[c(1:4, 8)], 
                         bal.stats$edital20.2016$balnoage[c(1:4, 8)]))

wald20 <- xtable(bal.stats$edital20.2016$wald)
bal.stats$edital20.2016$p.permutation #permutation p-value for Table C.5 edital 20

print(ed20,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c8.tex",sep=""))

print(wald20,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c5b.tex",sep=""))


################ Appendix D
#objects
load("Routputs/out-W1-attrition-joint.RData")
load("Routputs/out-W1-attrition.RData")

attritionW1 <- texreg(pol.attrition.stats, include.ci = F)

print(attritionW1,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d3.tex",sep=""))


#Joint hypotheses test
attrition_overall_covi <- pol.attrition.joint.stats[[1]] 
attrition_answered_covi <- pol.attrition.joint.stats[[2]]
attrition_overall_cov <- pol.attrition.joint.stats[[3]]
attrition_answered_cov <- pol.attrition.joint.stats[[4]]


joint_overall <- xtable(waldtest(attrition_overall_covi, attrition_overall_cov))
joint_within <- xtable(waldtest(attrition_answered_covi, attrition_answered_cov))

print(joint_overall,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d1a.tex",sep=""))

print(joint_within,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d1b.tex",sep=""))

########################## Appendix E

#Main ITT and CACE effects
load("Routputs/out-W1-estimates-pooled.RData")

#Getting estimates for Figure 1(a) and Table E1
allc <- cbind(out.all$out.itt.lm$out.lmc 
              ,out.all$out.itt.lin$out.linc)  
allc <- allc[1:10,c(1,2,3,6,5,7,8,9,12,11,10)] 
xallc <- xtable(allc[c(1:4),], digits=c(0,3,3,3,3,3,3,3,3,3,3,0))

#Table for Fig 1a
print(xallc,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-e1.tex",sep=""))

#Table for Fig 2
load("Routputs/out-W1-estimates-pooled-hetwait.RData")

allcwf <- out.waitenroll[[2]]
allcwf <- allcwf[,c(-4, -8)]
xallcw <- xtable(allcwf, digits=c(0,3,3,3,3,3,3,0,3)) 
label(xallcw) <- paste("tab:mod:waitenroll",sep="")
caption(xallcw) <- paste("Heterogenous Effects of Wait Since Enrollment (with controls)",sep="")
print(xallcw,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-e3.tex",sep=""))


######################### fig-e1a

the.labels.plot_tp <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                         "Attribution to\nLula/Dilma/PT", 
                         "Attribution to Paes/\nPezao/Temer/MDB") 


itts <- cbind(out.all$out.itt.lin$out.lin 
              ,out.all$out.itt.lin$out.linc)  
itts_tp <- itts[1:4,c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp
pdf(file = paste("Figures/fig-e1a.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(4,12,2,1))
plot(xs,c(-.75,-(nrow(tp)+2)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, 
     font = c(2, 1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1)) 
abline(v=0,lty=1,col=gray(0.7))
abline(h=c(ys)+0.5,lty=3,col=gray(0.7)) #adding two rows for comparison with W1
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT (via Lin)","ITT+Ctrls (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.015),xpd=NA)
dev.off()

############ Appendix F 

#Disaggregated plot
the.labels.plot_tp <-  c("Lula\nEvaluation", "Dilma\nEvaluation", "PT ID", 
                         "Paes\nEvaluation", "Pezão\nEvaluation", "Temer\nEvaluation")

itts <- cbind(out.all$out.itt.lm$out.lmc 
              ,out.all$out.itt.lin$out.linc)  
itts_tp <- itts[5:10,c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp
pdf(file = paste("Figures/fig-f1a.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(5,8,4,1))
plot(xs,c(-.75,-nrow(tp)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
     ,main="Recent Lotteries Survey")
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, 
     font = c(2, 1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1)) 
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT+Ctrls (via Fixed Effects)","ITT+Ctrls (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.0375),xpd=NA)
dev.off()


#Disaggregated plot without controls
the.labels.plot_tp <-  c("Lula\nEvaluation", "Dilma\nEvaluation", "PT ID", 
                         "Paes\nEvaluation", "Pezão\nEvaluation", "Temer\nEvaluation")

itts <- cbind(out.all$out.itt.lm$out.lm 
              ,out.all$out.itt.lin$out.lin)  
itts_tp <- itts[5:10,c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp
pdf(file = paste("Figures/fig-f2a.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(5,8,4,1))
plot(xs,c(-.75,-nrow(tp)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
     ,main="Recent Lotteries Survey")
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, 
     font = c(2, 1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1)) 
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT (via Fixed Effects)","ITT (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.0375),xpd=NA)
dev.off()


load("Routputs/out-W1-pred-wait.RData")

predictors <- pred_wait[[1]]
pred_int <-  pred_wait[[2]]

#Balance wait time
print(predictors,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-f1.tex",sep=""))

print(pred_int,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-f2.tex",sep=""))


############Appendix G

#Main ITT and CACE effects
load("Routputs/out-W1-estimates-pooled.RData")

the.labels.plot_tp <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                         "Attribution to\nLula/Dilma/PT", 
                         "Attribution to Paes/\nPezao/Temer/MDB")

caces <- out.all$out.cace.lm$out.lmcacec
caces <- caces[1:4,]
tp <- caces
pdf(file = paste("Figures/fig-g1a.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.8,.8)
par(mar=c(4,12,2,1))
plot(xs,c(-.75,-(nrow(tp)+2)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Complier Average Causal Effects"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp
     ,las=2,tick=F,cex.axis=1.1, font = 2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-6.5),col=gray(0.7))
abline(h=c(-5,-6)+0.5,lty=3,col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])
}
dev.off()


